
# This codes simulates positive semi-definite correlations matrizes of 
# nxn Dimension, based on the theorem 3.1 in
# Werner Huerlimann (2014): Cartesian and Polar Coordinates for the N-Dimensional Elliptope. Theoretical Mathematics & Applications 4(3):1-17.


huerlimann <- function(n, z){
  
  M <- matrix(NA, n, n)
  diag(M) <- 1
  for(i in 1:n-1){
    for(j in (i+1):n){
      M[i,j] <- runif(1, -1,1)   
    }
  }
  
  d <- n
  for(n in 1:d){
    for(i in 1:d){
      M[i,n] <- M[n,i]
    }
  }
  
  
  for(i in 1:(n-1)){
    
    M[i,n] <- z[i]
    
  }
  
  
  ### Generate empty nxn correlation matrix R
  
  R <- matrix(NA, n, n)
  diag(R) <- 1
  
  ### STEP 1 (Huerlimann 2014: 3.1)
  for(i in 1:(n-1)){
    
    R[i,n] <- M[i,n] 
    
  }
  
  ### STEP 2 (Huerlimann 2014: 3.2)
  for(i in 1:(n-2)){
    
    R[i,n-1] <- M[i,n] * M[n-1, n] + M[i, n-1]*sqrt((1- M[i, n]^2)*(1- M[n-1, n]^2))
    
  }
  
  ### STEP 3 (Huerlimann 2014: 3.3)
  for(k in 2:(n-2)){
    for(i in 1:(n-k-1)){
      
      
      p.1 <- M[i,n] * M[n-k, n]   
      
      prod.f  <- function(j){
        l <- (n-j+2):n
        prod(sqrt((1- M[i, l]^2)*(1- M[n-k, l]^2)), na.rm=T)  
      }
      prod.f <- Vectorize(prod.f, "j")
      
      j <- 2:k
      p.2 <- sum(M[i, n-j+1] * M[n-k, n-j+1] * prod.f(j) )
      
      l <- (n-k+1):n
      p.3 <- prod(sqrt((1- M[i, l]^2)*(1- M[n-k, l]^2)), na.rm=T) 
      p.3 <- M[i, n-k] * p.3
      
      
      R[i, n-k] <- p.1 + p.2 + p.3
      
    }
  }
  
  ##################################################################################################
  
  for(n in 1:d){
    for(i in 1:d){
      R[i,n] <- R[n,i]
    }
  }
  
  return(R)
  
}